

In [2]:

    
using Plots, ApproxFun, Interact
gr()  # Also works with plotly();

$$u'' + 2u = \cos(\omega t)\qquad u(0)=u'(0)=0$$

Slide to $\omega=\sqrt 2$ to see a resonance, where the solution grows



In [3]:

    
t=Fun(identity,[0.,1000.]) 
L=𝒟^2+2I  # our differential operator, 𝒟 is equivalent to Derivative()

@manipulate for ω=0.:.1sqrt(2):2sqrt(2)
    u=[ivp();L]\[0.;0.;cos(ω*t)]
    plot(u)
end









    











    



[Plots.jl] Initializing backend: gr






    Out[3]:

Bessel equation $$x^2 u'' + x u + (x^2 - \nu^2)u$$



In [4]:

    
x=Fun(identity,[1.,2000.])                                 
d=domain(x)
B=dirichlet()  
@manipulate for ν=0:.1:2000
	L=x^2*𝒟^2 + x*𝒟 + (x^2 - ν^2)   # our differential operator
   	u=[B;L]\[besselj(ν,first(d)),besselj(ν,last(d))]
    plot(u;ylims=(-.2,.2))
end









    











    Out[4]:

Exponentially ill-conditioned Lee & Greengard BVP $\epsilon u'' - x u' + u = 0$



In [5]:

    
x=Fun(identity)
B=dirichlet()

@manipulate for ϵ=0.0005:0.0005:0.4
    L=ϵ*𝒟^2 - x*𝒟 + I

     u=[B;L]\[1.,2.]
    plot(u;ylims=(-2.,3.))
end









    











    Out[5]:

Piecewise $$u'' + \omega u \begin{cases} 1 & |x| < c\cr 0 & {\rm otherwise}\end{cases} = 0,\qquad u(-1)=1,u(1)=0$$



In [6]:

    
@manipulate for ω=1:50000,c=0.1:0.01:1.
    f=Fun(x->abs(x)<c?1:0,[-1.,-c,c,1.])
    S=space(f)
    B=dirichlet(S)
    plot([B;𝒟^2+ω*f]\[1.,0.];ylims=(-0.5,1.1))
end

Forced Helmholtz equation $$u_{xx} + u_{yy} + k u = {\rm e}^{-(x-x_0)^2-(y-y_0)^2},\qquad \partial u = 0$$



In [11]:

    
plotly()  # switch to plotly to overcome bug in GR

d=Interval()^2
B=dirichlet(d)
Δ=Laplacian(d)

@manipulate for ny=10:200,x0=-1.:0.1:1.,y0=-1.:0.1:1.,k=-1000.:0.1:2000.
    f=Fun((x,y)->exp(-(x-x0)^2-(y-y0)^2),d)
    u=pdesolve([B;Δ+k*I],[zeros(∂(d));f],ny)
    contour(u)
end



In [12]:

    
d=Interval()^2
B=[dirichlet(d[1])⊗I;I⊗ldirichlet(d[2]);I⊗rneumann(d[2])]
Δ=lap(d)

@manipulate for ny=10:200,x0=-1.:0.1:1.,y0=-1.:0.1:1.,k=-1000.:1000.
    f=Fun((x,y)->exp(-(x-x0)^2-(y-y0)^2),d)
    u=pdesolve([B;Δ+k*I],[zeros(∂(d));f],ny)
    contour(u)
end

Dirichlet Helmholtz equation $$\Delta u + k u = 0, u(\pm 1,y)=u(x,\pm 1) = 1$$



In [13]:

    
d=Interval()^2
B=dirichlet(d)
Δ=lap(d)

@manipulate for k=-500.0:.001:2000.0,ny=10:200
    contour(pdesolve([B;Δ+k*I],ones(∂(d)),ny))
end

Convection diffusion $u_t = \epsilon u_{xx} + (B+C x) u_x, u(x,0)=e^{-20x^2}, u(-1,t)=u(1,t) = 0$



In [14]:

    
dx=Interval();dt=Interval(0,1.)
d=dx*dt
Dx=Derivative(d,[1,0]);Dt=Derivative(d,[0,1])
x=Fun(identity,dx)
@manipulate for ε=0.001:0.001:2.,B=-5.:0.1:5.,C=-5.:0.1:5.
    V=B+C*x    
    # Parentheses are a hack to get rank 2 PDE
    u=[timedirichlet(d);Dt-ε*Dx^2-V*Dx]\Fun(x->exp(-20x^2),dx)
    contour(u;xlims=(-1.,1.),ylims=(0.,1.))
end

Schrodinger equation $$i \epsilon u_t = -{\epsilon^2 \over 2} u_{xx} + x^2 u$$



In [15]:

    
dx=Interval(0.,1.);dt=Interval(0.0,0.54)
d=dx*dt

V=Fun(x->x^2,dx)

Dt=Derivative(d,[0,1]);Dx=Derivative(d,[1,0])

@manipulate for ϵ=0.005:0.002:0.3,ny=50:200
    u0=Fun(x->exp(-25*(x-.5)^2)*exp(-1.im/(5*ϵ)*log(2cosh(5*(x-.5)))),dx)
    L=1im*ϵ*Dt+.5*ϵ^2*Dx^2-V⊗1

    u=pdesolve([timedirichlet(d);L],u0,ny)
    contour(real(u))
end

Nonlinear BVP $$\epsilon u'' + 6(1-x^2)u' +u^2=1$$ $$u(-1)=u(1)=0$$



In [17]:

    
x=Fun()
u0=0.x
@manipulate for ε=0.001:0.001:1.,c=0.:0.1:1.
    N=u->[u(-1.)-c;u(1.);ε*u''+6*(1-x^2)*u'+u^2-1.]

    u=newton(N,u0)
    plot(u;ylims=(-1.,1.))
end